REPORT  DOCUMENTATION  PAGE 

Form  Approved 

0MB  NO.  0704-0188 

Public  reporting  Duraen  tor  this  colleciion  of  ntormaiion  is  estimated  to  average  1  hour  per  response,  mcluoing  the  time  tor  •"‘•'“'I®"® 

gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information  Send  comment  burden  «tm  *i2i^jSferson 

collection  of  information,  mduding  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  |®fvi«s^Dir^oi^  J«««rson 

Davis  Highway.  Suite  1204.  Arlington!^  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Project  {0704-0188).  Washington.  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

9/18/96  - 

4.  TITLE  AND  SUBTITLE 

Numerical  Solution  of  Nonlinear  Oscillatory  Multibody 

Dynamic  Systems 

5.  FUNDING  NUMBERS 

6.  AUTHOR{S) 

Jeng  Yen  and  Linda  Petzold 

7.  PERFORMING  ORGANIZATION  NAMES(S)  AND  ADDRESS(ES) 

University  of  Minnesota 

Department  of  Computer  Science 

4-192  EE/CS  Bldg,  200  Union  St  SE 

Minneapolis,  MN  55455 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSORING  /  MONITORING 

AGENCY  REPORT  NUMBER 

U.S.  Army  Research  Office 

P.O.  Box  12211 

Research  Triangle  Park,  NC  27709-2211 

11.  SUPPLEMENTARY  NOTES 

The  views,  opinions  and/or  findings  contained  in  this  report  are  those  of  the  author(s)  and  shoidd  not  be  construed  as 
an  official  Department  of  the  Army  position,  policy  or  decision,  unless  so  designated  by  other  documentation. 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  unlimited. 


1 9961 024  1 1 1 


13.  ABSTRACT  (Maximum  200  words) 

One  of  the  outstanding  problems  in  the  numerical  simulation  of  mechanical  systems  is  the 
development  of  efficient  methods  for  dealing  with  highly  oscillatory  systems.  These  type: 
of  sytems  arise  for  example  in  vehicle  simulation  in  modeling  the  suspension  system  or 
tires,  in  models  for  contact  and  impact,  in  flexible  body  simulation  from  vibrations  in 
the  structural  model,  and  in  molecular  dynamics.  Simulations  involving  high  frequency 
vibration  can  take  a  huge  number  of  time  steps,  often  as  a  consequence  of  oscillations 
which  are  not  physically  important.  The  components  causing  the  oscillations  cannot 
usually  be  eliminated  from  the  model  because  in  some  situations  they  are  critical  to  the 
simulation.  The  equations  of  motion  of  a  multibody  mechanical  system  are  described  by 
a  system  of  dif f erential-algrebraic  equations  (DAEs).  In  this  paper,  we  will  explore  two 
types  of  methods.  The  first  class  of  methods  damps  out  the  oscillation  via  highly  stable 
implicit  methods.  Even  in  this  relatively  simple  approach,  unforseen  problems  may  arise 
for  Newton  iteration  convergence,  due  to  the  nonlinearities.  The  second  class  of  methods 
involves  linearizing  the  system  around  the  smooth  solution.  The  linearized  system  can 
be  solved  rapidly  via  a  number  of  different  methods. 


14.  SUBJECT  TERMS 

constrained  dynamics,  multibody  systems,  differential-algebraic 
equations,  numerical  methods,  highly  oscillatory  systems. 


15.  NUMBER  IF  PAGES 
15 


16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  I  18.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT 
OR  REPORT  OF  THIS  PAGE  OF  ABSTRACT 

UNCLASSIFIED  UNCLASSIFIED  UNCLASSIFIED _ UL _ 


UNCLASSIFIED 


NSN  7540-01 *280-5500 


tff?0  QTJ.?'.LITY  MSPEGIS: 


_ UL _ 

standard  Form  298  (R«v.  2-89) 
PrMcrib^  by  ANSI  SlcL  239- 1 8 
298-102 


COMPUTER 


University  of  Minnesota 


INSTITUTE  OF  TECHNOLOGY 
UNIVERSITY  OF  MINNESOTA 

4-192  EE/CSci  Building 
200  Union  Street  SE 
Minneapolis,  Minnesota  55455 
612/625-4002 


TR  95-066 

Numerical  Solution  of  Nonlinear 
Oscillatory  Multibody  Dynamic 
Systems 

by  Jeng  Yen  and 
Linda  R.  Petzold 


Technical  Report 


Department  of  Computer  Science 
University  of  Minnesota 
4*192  EECS  Building 
200  Union  Street  SE 
Minneapolis,  MN  55455-0159  USA 


TR  95-066 

Numerical  Solution  of  Nonlinear 
Oscillatory  Multibody  Dynamic 
Systems 

by  Jeng  Yen  and 
Linda  R.  Petzold 


Numerical  Solution  of  Nonlinear 
Oscillatory  Multibody  Dynamic  Systems 


Jeng  Yen  * 
Linda  R.  Petzold  ^ 

August  24,  1995 


Abstract 

One  of  the  outstanding  problems  in  the  numerical  simulation  of  mechanical 
systems  is  the  development  of  efficient  methods  for  dealing  with  highly  oscilla¬ 
tory  systems.  These  types  of  systems  arise  for  example  in  vehicle  simulation  in 
modeling  the  suspension  system  or  tires,  in  models  for  contact  and  impact,  in 
flexible  body  simulation  from  vibrations  in  the  structural  model,  and  in  molec- 
ulcir  dynamics.  Simulations  involving  high  frequency  vibration  can  take  a  huge 
number  of  time  steps,  often  as  a  consequence  of  oscillations  which  aure  not  phys¬ 
ically  important.  The  components  causing  the  oscillations  cannot  usually  be 
eliminated  from  the  model  because  in  some  situations  they  are  critical  to  the 
simulation.  The  equations  of  motion  of  a  multibody  mechanical  system  are  de¬ 
scribed  by  a  system  of  differentiad-algebraic  equations  (DAEs).  In  this  paper, 
we  will  explore  two  types  of  methods.  The  first  class  of  methods  damps  out  the 
oscillation  via  highly  stable  implicit  methods.  Even  in  this  relatively  simple  ap¬ 
proach,  unforseen  problems  may  arise  for  Newton  iteration  convergence,  due  to 
the  nonlinearities.  The  second  class  of  methods  involves  linearizing  the  system 
around  the  smooth  solution.  The  linearized  system  can  be  solved  rapidly  via  a 
number  of  different  methods. 
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1  Introduction 


Much  recent  work  has  focused  on  the  development  of  numerical  methods  and  underly¬ 
ing  theory  for  the  solution  of  multibody  dynamic  systems  (MBS)  consisting  of  fast  and 
slow  subsystems  [17,  22].  These  types  of  systems  occur  frequently  as  initial  value  prob¬ 
lems  in  the  computer-aided  design  and  modeling  of  constrained  mechanical  systems, 
molecular  dynamics,  and  in  many  other  applications  [1,  31).  It  is  well-known  that 
the  characteristics  of  fast  or  slow  solution  are  determined  not  only  by  the  modeling 
aspects,  e.g.,  the  coefficients  of  stiffness  and  damping,  but  also  by  the  initial  conditions 
and  events  that  may  excite  stiff  components  in  the  system  during  the  simulation.  As 
an  example,  the  governing  equations  of  motion  of  a  mechanical  system  of  stiff  or  highly 
oscillatorv  force  devices  may  be  written  as  a  system  of  differential- algebraic  equations 
(DAE)  [5]: 

M{q)q  +  G'^{q)X-(f^  +  n  =  0  (la) 

g{q)  =0  (lb) 

where  q  =  [qi . qnf  is  the  generalized  coordinate,  g  =  ^  is  the  generalized  velocity, 

q-^  is  the  acceleration  and  A  =  [Ai . A,„]^  is  the  Lagrange  multiplier.  The  stiff 

or  oscillator}-  force  is  ft^  and  /"  includes  all  the  field  forces  and  the  external 

forces  which  are  non-stiff  compared  to  the  stiff  components,  e.g., 


'  ag  "  "  dq  ' 


(2) 


The  kinematic  constraints  are  g,  and  G  =  ||,  and  M  is  the  mass-inertia  matrix.  For 
the  stiff  force  components  in  (la),  we  assume  that 

//  =  -Biiq){KMq)  +  Gi^)  (3) 


where  77,  is  smooth  Vi  €  ,n/}  and  ,  and  Ki,  Ci  are  the  associated 

stiffness  and  damping  factors.  For  some  generalized  coordinate  sets,  the-  function  tj 
may  be  linear  or  even  identities,  e.g.,  for  instance  rji  =  g^,  for  some  i  €  {1, ....  n/}  and 
ii  €  {!,... ,n}.  Wffien  the  components  of  the  coefficient  matrices  Ki  and  C,  become 
large,  these  force  components  may  cause  rapid  decay  or  high  frequency  oscillation  in 
the  solution  of  (1).  The  purpose  of  this  article  is  to  study  these  systems  and  their 
numerical  solution.  In  this  notation,  the  stiff  force  term  in  (la)  can  be  written  as 

P  =-B{q){Kv[q)  +  CB[q)v).  (4) 


To  demonstrate  the  problem  of  oscillation  and  the  recent  developments  in  this  area, 
we  present  two  examples:  a  stiff  pendulum  amd  a  2D  bushing  problem.  The  former 
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is  a  ver>-  simple  example  of  a  type  of  system  often  seen  in  modeling  molecular  dy¬ 
namic  systems,  and  the  latter  is  a  general  form  of  modeling  force  devices  in  multibody 
mechanical  systems. 

Stiff  pendulum 

In  Cartesian  coordinates,  a  simple  stiff  pendulum  model,  with  unit  mass  and  gravity, 
may  be  expressed  as 


0 

=  X  —  u 

(5a) 

0 

=  y-v 

(5b) 

0 

=  u  +  xA 

(5c) 

0 

=  i’  +  yA  —  1.0 

(5d) 

y/x'^  +  y'^  —  1 .0 
y/x'^  +  y^ 

(5e) 

where  the  stiff  spring  of  natural  length  1.0  and  stiffness  ji  is  attached  to  the  center 
of  mass  of  the  pendulum.  Preloading  the  spring  by  using  c  =  \/10"^,  the  initial 
condition  (xo,yo)  =  (0.9. 0.1)  and  the  zero  initial  velocity  {uq,  t’o)  =  (0,0),  the  results 
of  the  states  (z,y,u,t’)  in  the  0  to  10  second  simulation  are  shown  in  Fig.  1.  The 
corresponding  eigenvalues  of  the  underlying  ODE  of  (5),  i.e.,  substituting  (5e)  into 
(5c.  5d),  are  illustrated  in  Fig.  2,  where  the  3D  figures  contain  all  the  eigenvalues  on 
the  complex  plane  drawn  along  the  time-axis.  The  dominant  pair  of  eigenvalues  in  the 
example  are  as  shown  in  Fig.  2.  As  e  0,  the  pair  of  eigenvalues  approaches 
±00  along  the  imaginary  axis.  The  other  pair  of  eigenvalues  oscillates  on  the  complex 
plane  with  the  amplitude  and  frequency  approaching  ±oo.  Decreasing  e  to  \/l0"=,  the 
eigenvalues  of  the  underlying  ODE  of  (5)  are  times  the  magnitude  of  those  in  Fig.  2, 
and  the  oscillating  pair  increases  its  frequency  proportional  to  the  size  of  e. 


s — ; — t 


\ — ; — t — ; — i — s — J — i — s — 

l* — 

\ — ; — ; — ! — .  — i — s — ; — » 

Figure  1:  Stiff  Pendulum  in  Cartesian  Coordinates 

Lubich  [17]  shows  that  the  numerical  solution  by  a  class  of  Runge-Kutta  methods 
of  stiff  mechanical  systems  of  a  strong  potential  energy,  e.g.,  stiff  spring  force  such  as 
the  stiff  pendulum  (5),  converges  to  the  slowly  varying  part  of  the  solution,  with  the 
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Figure  2:  Eigenvalues  of  Stiff  Pendulum  in  Cartesian  Coordinates,  epsilon  -  lOe-1.5 

stepsize  independent  of  the  parameter  e  in  (5).  Reich  [22]  extends  the  principle  of 
slow  mani fold  [9,  15]  to  DAE  of  MBS  with  highly  oscillatory  force  terms.  Algebraic 
constraints  corresponding  to  the  slow  motion  were  introduced  with  a  relaxation  pa¬ 
rameter  to  preserve  the  slow  solution  while  adding  flexibility  to  it  in  the  slow  manifold 

approach. 

It  is  not  clear  that  a  slow  solution  appears  in  the  above  example.  In  fact,  we  can 
only  identify  the  slow  solution  of  (5)  using  a  proper  nonlinear  coordinate  transforma¬ 
tion.  In  polar  coordinates  (r,^),  we  obtain  the  equations  of  motion  of  (o): 


0  =  r-2 

0  =  ^-u; 

0  =  i-t-rw^-l- ^(r- 1) -sin^  (6c) 

0  =  w --{2zu)  -  cos9)  (6d) 

r 


where  {z,w)  is  the  velocity.  In  the  0  to  10  second  simulation,  using  the  same  initial 

conditions  as  the  previous  examples,  i.e.,  r©  =  y/il  +  vh  =  arctan  and  (zq,uio)  = 
(0,0),  where  (xci/o)  is  the  initial  position  of  (5),  we  obtain  the  solution  in  Fig.  3.  It 
is  clear  that  the  solutions  of  (r,r)  represent  the  fast  motion,  and  the  solutions  of 
(9,  w)  are  the  slow  ones.  The  eigenvalues  along  the  solution  trajectory  are  presented 
in  Fig.  4.  Note  the  dominate  eigenvalues  are  of  the  same  magnitude  as  those  in  (5), 
see  Fig.  2.  This  is  because  the  coordinate  transformation,  x  =  rcosd,  y  =  rsind,  is 
linear  with  respect  to  the  fast  moving  r.  The  eigenvalues  of  (6)  with  c  =  Vl^  are 
similar  to  the  comparison  in  the  Cartesian  formulation,  i.e.,  we  obtain  the  eigenvalues 
of  10  times  magnification.  However,  the  eigenvalues  corresponding  to  the  slow  motion 
have  near  zero  imaginary  parts,  therefore,  the  oscillations  along  the  imaginary  axis  of 
eigenvalues  3,4  in  Fig.  4  remain  insignificant. 

Although  there  are  ongoing  developments  to  extend  the  results  of  Lubich  to  multi¬ 
stage  multistep  methods  [23],  and  impressive  application  of  the  slow  manifold  tech¬ 
nique  in  some  molecular  dynamic  models,  it  is  not  clear  that  these  results  may  apply 
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Figure  3;  Stiff  Pendulum  in  Polax  Coordinates 


Figure  4:  Eigenvalues  of  Stiff  Pendulum  in  Polar  Coordinates,  epsilon  =  lOe-1.5 

directly  to  all  the  types  of  oscillatory  components  in  MBS.  As  indicated  in  [17],  the 
representation  of  stiff  or  oscillatory  components  in  an  appropriate  coordinate  system 
of  MBS  is  not  always  possible,  i.e.,  the  constraiints  associated  with  the  stiff  or  oscil¬ 
latory  potential  force  can  be  difficult  to  obtain  in  general.  Nevertheless.  for  (1),  an 
approximation  of  the  dynamics  of  such  local  coordinates  can  be  obtained  for  the  os¬ 
cillatory  components  in  the  form  of  (4).  We  are  also  concerned  with  the  convergence 
of  New'ton’s  method  for  the  numerical  solution  with  large  stepsize  of  (1),  which  may 
be  an  obstacle  in  obtaining  efficient  numerical  solution  of  an  oscillatory  MBS  in  either 
of  the  above-mentioned  approaches. 

Bushing  force 

We  have  been  studying  more  general  MBS  of  nonlinear  oscillatory  components  such  as 
a  bushing  force,  which  is  often  used  in  modeling  vehicle  suspension  systems.  Different 
from  the  linear  spring,  this  element  is  usually  an  anisotropic  force,  i.e.,  it  has  different 
spring  coefficients  along  the  principle  axes  of  the  bushing  local  coordinate  frame.  The 
bushing  force  between  body-i  and  body-j  may  be  defined  using  the  relative  displacement 
d,j,  its  time  derivative  d,j,  and  the  relative  angle  Oij  and  its  time  derivative  6ij  of  two 
body-fixed  local  coordinate  frames  at  the  bushing  location  on  two  bodies.  Using  the 
vectors  s'  and  s'  representing  the  bushing  location  in  the  body-ts  and  body-j's  centroid 
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local  coordinate  systems,  respectively,  we  have 


where  the  orientation  transformation  matrices  Ai  and  Aj  are 


j  ^  _  f  cosd,  -sind, 

.4,  =  cos^,- 


A  A/a  \  COS(7j  —  smc7j 

Aj  =  A(S,)  =  [  J 

and  [x,,  t/i,  5,]  and  [xj,  yj,  dj]  are  Cartesian  coordinates  at  body-fixed  frames.  The 
bushing  force  fb  can  then  be  written  as 


A  =  f  1  =  A 


1  -I-  4  f  ^  ^1  4^d 

Q  J^y  Ai  dij  +  Ai  Q  Ai  C,j 


and  the  applied  torque  is 

Tb  =  k^dij  +  ^ij  (^) 

where  ju'ij  =  ^,  A:^,  and  are  the  spring  coefficients  associated  with  the  x,  y, 
and  9  coordinates,  and  c®,  c*',  and  c®  are  the  corresponding  damping  coefficients. 


A  simple  example  may  be  obtained  from  this  model  using  unit  mass-inertia  and 
gravity,  grounding  the  first  body,  and  setting  the  bushing  location  on  the  second  body 
to  5'  =  [-5,  0].  A  bushing  element  with  no  damping  attached  at  the  global  position  of 
[5,  0]  yields 

0  =  x-k^{^-x  +  ^^)  (10a) 

0  =  y  +  kHy-^^)  +  l  (10b) 

0  =  d  +  k^9-—k^{--x  +  —)-—ky{y-—).  (10c) 

It  is  easy  to  see  from  (10)  that  the  local  eigenstructure  of  the  system  may  change 
rapidly,  depending  on  the  size  of  the  stiffness  coefficients. 

Using  the  initial  values  of  (x,  y,  9)  =  (1.1, 0.1, 0.0)  with  {k^,  fc*',  k^)  =  (10'’,  10^,  10^), 
the  solution  of  (10)  exhibits  high  frequency  oscillation  for  all  coordinates,  as  shown  in 
Fig.  5.  Solving  the  eigenvalue  problem  of  (10)  at  each  time  step  yields  three  pairs  as 
illustrated  in  Fig.  6. 

Many  methods  for  efficient  solution  of  oscillatory  dynamic  systems  are  predicated 
on  a  nearly  linear  form  of  the  equation.  For  example,  the  method  of  averaging  [4] 
requires  the  linear  part  of  the  oscillation  equations  of  motion  to  be  dominant,  and  the 
mode-acceleration  method  for  structural  dynamics,  which  eliminates  higher  modes 
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Figure  5:  Bushing  Problem  in  Cartesian  Coordinates 


Figure  6:  Eigenvalues  of  Bushing  Problem 

in  the  computation  of  the  mode-displacement  solution  [29],  is  based  on  the  time- 
invariant  eigenvalues  of  the  structural  dynamic  equations.  Our  aim  is  to  treat  the 
class  of  general  nonlinear  stiff  and  oscillatory  forces  represented  in  the  MBS  of  (1). 
One  approach  is  based  on  the  study  of  a  class  of  MBS  DAE  solvers  [27]  and  the  ener^ 
dissipative  method,  which  may  damp  out  the  oscillation  where  the  amplitude  is  small. 
The  other  approach  is  to  approximate  the  system  by  a  nearly  linear  system,  whose 
solution  approximates  the  solution  to  the  original  system. 

The  main  technical  problem  in  either  approach  arises  from  the  fact  that  for  these 
nonlinear  oscillating  problems,  the  eigenstructure  of  the  local  Jacobian  varies  rapidly, 
on  the  scale  of  the  high-frequency  oscillation,  in  a  neighborhood  of  the  smooth  solution. 
We  will  used  the  fact  that  the  eigenstructure  varies  slowly  along  the  smooth  solution 
to  construct  efficient  numerical  methods  and  accurate  approximate  solutions. 


2  Damping  the  oscillation 

Given  the  possibility  of  a  rapidly  changing  local  eigenvalue  structure,  perhaps  the 
simplest  strategy  is  to  consider  damping  the  oscillation  whenever  it  is  not  important  via 
highly  stable  implicit  numerical  methods.  Since  the  amount  of  damping  is  controlled 
by  the  time-step,  and  automatic  stepsize  selection  increases  the  time-step  whenever  the 
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solution  is  slowly- vamng  (i.e.  if  the  amplitude  of  the  oscillation  is  small  m  comparison 
with  the  local  error  tolerances),  the  stepsize  is  increased  when  the  oscillation  is  no 

longer  important. 

In  recent  work  [27]  we  have  considered  the  solution  of  mechanical  systems  with 
high  frequency  vibrations  via  this  type  of  technique.  In  experiments  wuth  the  bushing 
problem  (10)  solved  directly  by  low-order  BDF  methods,  we  found  that  the  methods 
experienced  severe  problems  with  Newton  convergence.  To  overcome  these  problems, 
we  proposed  a  coordinate-split  (CS)  formulation  of  the  equations  of  motion,  and  a 
Newton-type  iteration  for  solving  the  coordinate-split  equations  at  each  time  step. 
The  coordinate-split  formulation  eliminates  problems  due  to  obtaining  an  accurate 
predictor  for  the  Lagrange  multiplier  variables  because  these  vanables  are  no  longer 
present  in  the  computation.  We  found  that  the  coordinate-split  formulation  worked 
well  for  several  test  problems  involving  mechanical  systems  with  high  frequency  os¬ 
cillations.  However,  for  problems  with  very  high-frequency  oscillations,  there  are  still 
difficulties  in  Newton  convergence  with  this  method. 

The  Jacobian  matrix  for  solving  the  nonlinear  equations  of  the  coordinate-split 
formulation  at  each  time  step  involves  several  terms  which  are  complicated  to  com¬ 
pute  and  which  are  small  at  the  solution  of  the  nonlinear  system.  These  are  terms  of 
second-order  which  correspond  to  the  derivative  of  the  projection  operator  onto  the 
constraints.  Awav  from  the  smooth  solution,  these  terms  are  highly  oscillatory.  By 
neglecting  these  terms,  we  found  that  the  resulting  New-ton-type  method  converged 
much  faster  for  oscillating  test  problems  like  the  bushing  problem.  We  called  the  re¬ 
sulting  method  the  modified  coordinate- split,  or  CM  method.  In  [28],  convergence  of 
the  CM  iteration  is  analyzed,  and  the  improved  convergence  for  oscillatory  multibody 
svstems  is  explained.  Intuitively,  by  neglecting  these  terms  the  CM-iteration  yprox- 
imates  the  Jacobian  along  the  smooth  solution,  thus  yielding  more  reliable  Newton 

directions. 

The  modified  coordinate-split  (CM)  method  performed  extremely  well  in  numer¬ 
ical  experiments  described  in  [27].  The  constraints  g{p)  =  [gu-.gs]  of  a  two-body 
pendulum  may  be  wTitten  as 


gi  = 

(11a) 

92  =  yi 

(11b) 

C 

93  =  ^1 

(11c) 

p4  =  (xi  -  X2)^  +  (yi  ”  92)^  -  1 

(lid) 

♦ 

95  =  ^2 

(lie) 

where  x,  and  j/,-,  i  =  1,2  are  Cartesian  coordinates  of  the  center  of  mass  of  body  i, 
and  d,  is 'the  orientation  coordinate  of  the  body  centroid  reference  coordinate^systein, 
and  tile  length  of  the  pendulum  is  1.  Applying  the  bushing  force  (8)  with  [fc^,P,fc  ] 
=  [1000, 1000, 1000]  and  [c*,  c^,  c®]  =  [10. 10, 10]  to  the  pendulum,  small  oscillations  of 
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the,  numerical  solution  appear.  Using  the  initial  values  q 

and  V  =  [0,0,0,-6.75e-5,-4.5444e-3],  numerical  results  from  the  BDF  code  DAbbL  [21] 
are  contained  in  Table  U  in  which  error  test  failures  {etf  -  -s)  and  convergence  test 
failures  (ctf  -  s)  are  listed.  We  denoted  by  CS  the  coordinate-splitting  formulation. 
LG  the  stabilized  index-2  formulation  proposed  by  Gear  [10],  CM  the  coordinate- 
split  form  using  a  modified  iteration  matrix  with  the  second-order  derivative  terms 
omitted,  and  LM  the  modified  LG  using  the  new  predictor  of  the  multipliere  by  the 
CS  method.  Using  simplified  Newton  iterations  and  the  corresponding  modified  local 
error  estimate,  CS.  CM.  LG  and  LM  obtain  consistent  results. 


Method 

TUT' 

step 

etf  -  s 

ctf  -  .s 

CS 

~W^ 

62 

156 

48 

0 

13 

CM 

10'® 

62 

156 

48 

0 

13 

LM 

10'® 

62 

154 

48 

0 

13 

LG 

10'® 

59 

141 

46 

0 

12 

CS 

lO''' 

77 

193 

65 

1 

16 

CM 

lO"* 

77 

193 

65 

1 

16 

LG 

10'^ 

61 

136 

27 

1 

5 

LM 

10'^ 

77 

193 

65 

1 

16 

CS 

10'® 

87  1 

215 

54 

0 

13 

CM 

10'® 

87 

215 

54 

0 

13 

LG 

10'® 

108 

259 

77 

1 

21 

LM 

10'® 

87 

215 

54 

0 

13 

CS 

10'® 

138 

343 

97 

1 

25 

CM 

10'® 

138 

343 

97 

1 

25 

LG 

10'® 

131 

308 

65 

0 

16 

LM 

.  10'® 

138 

343 

97 

1 

25 

Table  1;  Simple  Pendulum  with  a  Bushing  Force.  Spring  Constant  -  10^ 


To  see  the  effect  of  more  severe  oscillation,  we  increased  the  spring  constant  of 
the  bushing  to  10®.  Time  steps  of  these  methods  selected  by  DASSL  are  shown  in 
Figure  7.  Clearlv,  CM  took  much  larger  steps  than  the  other  methods.  Moreover,  if 
the  spring  constant  is  increased  to  10®,  we  found  severe  convergence  problems  for  LG 
CS  and  IM;  the  results  are  contained  in  Table  2.  Further  details  on  the  numerical 
experiments  are  given  in  [27]. 


3  Smooth  Linearization 

Often  in  multibody  systems  the  components  exhibiting  high  frequency  oscillation  result 
from  the  potential  forces  induced  by  material  deformations.  In  flexible  multibody 
dynamic  systems,  for  example,  there  are  usually  nonlinear  transformations  applying 
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Method 

TUT 

time 

step 

/  -  s 

j  -  s 

etf  -  s 

ctf  -  s 

C5,„ 

io-'‘ 

0-0.1 

2252 

4901 

3361 

1 

1120 

CM,„ 

lO-'* 

0-0.1 

20 

40 

7 

0 

0 

LG,n 

10-^ 

0-0.1 

5267 

10536 

7899 

0 

2633 

LMsn 

lO-'* 

0-0.1 

2251 

4900 

3360 

1 

1120 

Table  2:  Results  of  Bushing  Problem,  Spring  Constant  =  10®,  Damping  =  10 


Spring  »  10«5  0*fnpr»g  *  10,  TOL  * 

O.OSr  CM 


Figure  7;  Time  Steps  Used  in  Solving  the  Bushing  Problem,  Spring  Constant  =  10® 
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to  the  internal  and  inertia  forces  of  the  flexible  structure  to  obtain  equations  of  motion 
in  the  generalized  coordinate  space  [29] .  The  resulting  system  of  equations  contains 
nonlinear  high  frequency  oscillatory  forces. 

One  approach  to  solving  the  high  frequency  oscillation  problem  is  to  carry  out 
modal  analysis  and  then  eliminate  the  higher  modes,  since  lower  modes  may  preserve 
the  slowly  vairying  part  of  the  solution  [8, 13).  For  example,  the  extreme  high  modes  of  a 
structure  are  often  rejected  in  modeling  flexible  effects  of  mechanisms,  since  the  details 
of  the  oscillating  solution  are  not  so  important  as  the  long-term  solution  behavior.  A 
similar  approach  has  been  developed  in  recent  work  on  molecular  dynamics  simulation 
[31].  However,  due  to  nonlinear  oscillatory  forces  in  the  multibody  formulation,  the 
modal  analysis  needs  to  be  carried  out  at  each  time  step  to  resolve  the  rapidly  varying 
local  eigenstructure  of  the  system,  resulting  in  very  costly  computations. 

Another  approach  is  to  resolve  the  oscillation  efficiently  via  look-up  tables.  This 
idea  is  frequently  used  in  a  real-time  simulation  environment.  One  such  example  is  the 
modeling  of  contact  compliance  in  rigid  body  simulation,  where  the  localized  oscillation 
may  be  of  interest.  Applying  linear  constitutive  laws  to  modeling  the  contact  com¬ 
pliance,  e.g.,  the  elastic  half  space  theory,  Boussinesq’s  influence  functions  and  Hertz’ 
contact  model  leads  to  linear  spring  forces  between  contact  bodies  [26,  14,  18].  The 
spring  coefficients  may  be  very  Itirge  since  the  contact  deformations  are  small  compared 
to  the  gross  motion  of  the  contacting  bodies.  The  advantage  of  using  table-look-up 
is  efficiency,  however  it  is  not  clear  how  the  variable  stepsize  and  order  numerical 
integration  should  interact  with  the  tables  to  maintain  efficiency  and  accuracy. 

In  a  numerical  method  such  as  multistep  or  Runge-Kutta,  which  are  based  on  ap¬ 
proximating  the  solution  locally,  the  stepsize  must  be  chosen  very  small  to  resolve  the 
high-frequency  oscillation  in  the  system.  Moreover,  due  to  the  nonlinear  transforma¬ 
tion  that  places  oscillating  components  in  the  space  of  the  generalized  coordinates, 
e.g.,  in  the  form  of  (1),  the  numerical  method  may  become  ineffective  since  the  eigen¬ 
structure  may  change  rapidly  as  shown  in  the  previous  examples.  Our  goal  of  treating 
large-scale  MBS  w’ith  highly  oscillatory  components  (1)  is  to  develop  numerical  meth¬ 
ods  that  approximate  the  high  frequency  components  properly. 

Modal  analysis  in  structural  dynamics  is  well-developed  and  implemented  in  pro¬ 
duction  software  [8,  19].  As  shown  in  [7],  combining  structural  dynamic  subsystems 
with  the  DAE  of  MBS,  high  frequency  nonlinear  oscillating  solutions  may  occur.  Based 
on  the  solution  of  this  class  of  nonlinear  oscillating  problems,  we  propose  a  new  ap¬ 
proach  which  is  based  on  linearizing  the  oscillating  components  around  the  smooth 
solution. 

Under  the  assumption  that  the  oscillatory  forces  are  of  the  form  (3),  we  can  choose 
the  local  coordinates  q  =  [gi, ...,  9n/]^  and  the  velocity  v  =  §  such  that 

Q-^liq)  =  0  (12a) 
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V  -  B{q)^v  =  0 


(12b) 


at  each  force  component.  Note  that  the  corresponding  potential  energy  becomes 

=  fKq  ■  (13) 

where  |1/C||  is  the  dominant  term  of  the  highly  oscillatory  forces.  The  oscillations  in 
the  solution  of  q  are  coupled  with  those  of  the  local  coordinate  q.  Since  the  potential 
(13)  is  quadratic  in  q,  the  dynamic  equations  of  q  contain  a  (nearly)  linear  term  that 
characterize  the  local  oscillations  due  to  the  forces. 

To  derive  the  equations  of  motion  for  q,  we  can  eliminate  A  in  (1)  using  the 
acceleration  equations, 

G{q)q  = 

obtained  from  twice  differentiating  (lb).  Let  M{q)  be  nonsingular  and  =  B{q)Kq, 
solving  for  q  from  (la)  and  (14)  yields 

q  ^  M-‘(/  -  G’'(GM-‘G’')-*GM-‘)(/”  +  n  +  (15) 


Differentiating  (12b)  with  respect  to  time  yields 


^T■.  -rdB^  ■ 

A'*  >-  — 


q-B{qrq-e-^Q  =  0. 


(16) 


(17) 

(18) 


Substituting  (15)  into  (16)  we  obtain 

'^  +  Kq  =  f{q/q,t) 

in  which  we  denote  M~^K  =  K  such  that 

M-l  =  B{qfMq)B{q) 

where  A(g)  =  -  G'^{GM-^C^)-^GM-^).  According  to  the  assumption  (2), 

the  right-hand  side  of  (17) 

/  =  7®-5^A5/"  (19) 

where 

7^  =  G^GM-^G^r^'r  + 

satisfies  - 

W^W  «  (20) 

Thus,  the  solution  of  (17)  has  an  asymptotic  expansion  with  a  dominant  first  term 
for  a  fixed  q  [4],  i.e.,  a  smooth  solution  of  q  can  be  expressed  as 


9  = 

i=o 


(21) 
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where  Qj  are  smooth  functions  and  He/CH  «  1. 

As  shown  in  the  stiff  pendulum  and  the  bushing  examples,  the  nonlinear  oscillatory, 
force  in  (1)  can  be  written,  using  proper  local  coordinates,  in  a  neairly  linear  form.  In 
the  stiff  pendulum,  the  proper  choice  of  the  local  coordinate  is  the  polar  system.  In 
the  bushing  example,  this  is  achieved  using  the  local  displacement  and  relative  angle 
at  the  bushing  reference  frame.  Using  this  approach,  the  general  form  of  flexible  MBS 
can  be  written  as 


=  0  (22a) 

z'  +  H{y)z  +  h{y,z,t)  =  0  (22b) 

where  y  =  [g,g]  and  2  =  [5,^],  and  (22a)  is  often  a  system  of  differential-algebraic 
equations.  Note  that  (22)  usually  contains  a  large  number  of  equations,  i.e.,  there  are 
many  of  the  free  vibration  modes. 

Our  objective  is  to  develop  a  method  that  takes  large  time-steps  relative  to  the 
high-frequency  oscillations.  Several  methods  are  currently  under  investigation.  One 
approach  is  to  solve  a  linearization  of  (22)  directly  by  damped  numerical  niethods.  The 
linearization  is  carried  out  aroimd  a  smooth  solution  y  such  that  y  =  [g,  g]  is  near  the 
equilibrium  of  the  highly  oscillatory  forces,  i.e.,  77(g)  ss  0  and  B^{q)q  «  0.  To  obtain 
the  numerical  solution,  we  apply  the  above-mentioned  CM  method  to  (22a),  and  the 
Lanczos  types  or  Amoldi  adgorithm  to  the  resulting  linear  system.  These  methods 
for  nonsymmetric  linear  systems  with  proper  modifications,  for  example  in  [12],  are 
used  to  acquire  a  stable  reduced  order  model  for  the  linearized  form  of  (22).  Another 
possibility  is  to  apply  modal  analysis  to  (22b)  along  the  smooth  solutions,  yielding  a 
reduced  system  of  (22b)  that  is  described  by  the  Ritz  vectors  or  the  Lanczos  vectors 
[25.  16,  30,  24,  19].  The  key  idea  is  to  obtain  a  smooth  linearization  by  choosing  the 
proper  local  coordinates  and  linearizing  along  the  smooth  solution. 
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